Turing patterns on networks 
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Turing patterns formed by activator-inhibitor systems on networks are considered. The linear 
stability analysis shows that the Turing instability generally occurs when the inhibitor diffuses 
sufficiently faster than the activator. Numerical simulations, using a prey-predator model on a scale- 
free random network, demonstrate that the final, asymptotically reached Turing patterns can be 
largely different from the critical modes at the onset of instability, and multistability and hysteresis 
are typically observed. An approximate mean- field theory of nonlinear Turing patterns on the 
networks is constructed. 
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PACS numbers: 82.40.Ck, 89.75.Fb, 87.23.Cc 

Turing instability in activator-inhibitor systems is one 
of the most important classical concepts of nonequilib- 
rium pattern formation [1], with diffusion destabilizing a 
uniform stationary state of the system and leading to for- 
mation of stationary spatial patterns. Turing patterns in 
ordinary continuous media have been extensively inves- 
tigated in theoretical [2] and experimental 0, 0] studies. 
Their properties have also been studied for spatially het- 
erogeneous systems [5[ and in the presence of fractional 
diffusion [6]. 

In this Letter, we analyze Turing patterns exhibited 
by activator-inhibitor systems on networks. In such sys- 
tems, the nodes are populated by activator and inhibitor 
species which are diffusively transported over the edges 
that form a network. As examples, multicellular sys- 
tems [7], networks of coupled chemical reactors [8], and 
transportation networks that facilitate spreading of an- 
imals or insects can be mentioned. If several 
species are present, their mobilities on a network may 
significantly differ, and a situation is possible where the 
inhibitor species is much more easily transported over 
a network than the activator. Previously, it has been 
shown by linear stability analysis that the uniform sta- 
tionary state in such networks can be unstable 
However, detailed numerical and analytical investigations 
have so far been performed only for lattices [7[ and for 
small network systems [8|- 

Here, we provide a complete linear stability analy- 
sis valid for any activator-inhibitor system and perform 
a detailed numerical study of Turing patterns in the 
Mimur a- Murray model of interacting prey-predator pop- 
ulations pjj . Our numerical results reveal that the actual 
nonlinear behavior of network-based activator-inhibitor 
systems may show significant differences from the pre- 
dictions of the linear stability analysis, with coexistence 
of many different structures and hysteresis being typi- 
cally observed. As we further show, the properties of 
network Turing patterns can often be well understood in 
the framework of a mean-field approximation. 

Activator-inhibitor systems on networks are described 



by equations 



Ui(t) =f{Ui, V z ) + e(\7 2 U) u 
V i (t)=g(U h V i ) + ae{V 2 V) h 



(1) 



where Ui and Vi are concentrations of the activator and 
the inhibitor species on node i = 1, • • • ,iV, and func- 
tions f(U,V) and g(U, V) determine their local dynam- 
ics on a single node. In these equations, V 2 represents 
the diffusion (or Laplacian) operator on the network de- 
fined as (V 2 ^); = E 



N 



i Aij(Uj 



Ui) and (V 2 V), = 

J2f=i Aij(Vj — Vi), where (Aij) is the adjacency matrix 
of the network. All nodes and edges in the network 
have the same properties. The parameter e is the dif- 
fusion constant of the activator, and a is the ratio of 
the inhibitor diffusion constant to that of the activator. 
Note that the action of the network diffusion operator 
can also be expressed as (V 2 U)i = Y^j=i^ij^j where 
(Lij) is the Laplacian matrix of the network, defined as 



L 



ki5ij with k{ 



E^Li Aij representing the 



ij — 

degree of the node i. 

In absence of diffusion (e = 0), each network element 
has a linearly stable stationary state (U(°\V^), which 
satisfies f(U(°\vW) = and g(U^\V^) = 0. Since U 
is the activator and V is the inhibitor, we assume f u = 
9f/dU\ iu(0)jV( o) ) > 0, f v = df/dV\ iu(0)y( o) ) < 0, g u = 
9g/dU\ {u{Q)y{Q)) > 0, and g v = dg/dV\ {u(0) y (0)) < 0. 
When diffusion is turned on (e > 0), the uniform solution 
{U u Vi) = (U®\V<®) still satisfies Eqs. Q. However, 
this uniform solution may now become unstable due to 
the difference in diffusional mobilities of the activator and 
inhibitor species. 

The linear stability analysis is naturally performed 
by using eigenvalues A^ and eigenvectors = 

(V 2 0(«)), 



N ) of the Laplacian matrix, satisfying 



1, 



,N). 



Since is a real symmetric matrix, all its eigenvalues 
are real, and its eigenvectors are mutually orthogonal. 
We sort indices {a} in the decreasing order of the eigen- 
values, so that = A^ > A( 2 ) > • • • > A W holds. 
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FIG. 1: (color online) Linear growth rates in the Mimura- 
Murray model on the scale-free network, with Re A^ plotted 
against — eA^ for e — 0.12 and three different values of the 
ratio a of the diffusion constants. Crosses show actual discrete 
values of the growth rates for a = 15.6. 



Linearized equations describing evolution of small per- 
turbations (iii,Vi) to the uniform state (U^°\V^) are 
given by 



Ui(t) 
Vi{t) 



f u Ui + f v Vi + e(V 2 u)i, 



g u Ui 



JvVi 



(2) 



By expanding (iii,Vi) over the Laplacian eigenvectors, 

(ui,Vi) = ELi(^ (a) ^ M ) ex P [^ a) t]4 a) where u(a) 
and are the expansion coefficients, Eqs. (j2j) are trans- 
formed into N independent linear equations for different 
modes (u^ a \v^). The linear growth rate A^ of the a- 
th mode is determined from the characteristic equation 
{X (a) -fu- eA(°) }{ A(°> -g v - aeA^ } - f v g u = 0, which 
yields 

A (a) = (/„+5„ + (l + <7) e A( a ) 

±\/4/i + (fu ~gv + (l- a)eA(«)) 2 ^) / 2. (3) 

When Re A^ is positive, the a-th mode is unstable. 
The instability takes place when one of the modes begins 
to grow, namely, when Re A^ = for some a = ciq 
and Re A^ < for all other modes. For the Turing 
instability, the condition Im A^ =0 should also hold. 

Graphically, all growth rates A^ lie on the curve 
A = F(A), which is determined by Eq. ^ (see Fig. [T]). 
When the ratio a of the diffusion constants is var- 
ied, this curve touches the horizontal axis at <j c = 

{fu9v ~ 2fv9u + 2y/fv9u (fv9u ~ f u g v )^j / fl, and some 
of the modes may therefore become unstable for a > <j c . 
For this to occur, the eigenvalues of the discrete Lapla- 
cian spectrum of the considered network must be actually 
present near the maximum of the curve. Therefore, the 
condition a = a c provides only the lower bound for the 
Turing instability boundary. 



Thus, the condition of the Turing instability on net- 
works is similar to the respective condition for the con- 
tinuous media. The role of plane waves is now played by 
Laplacian eigenvectors of the considered network, and, 
instead of wavenumbers, the respective Laplacian eigen- 
values are important. Note that the above linear stabil- 
ity analysis is general; it holds for any activator-inhibitor 
model and for any network architecture. The critical 
mode always corresponds to a certain Laplacian eigen- 
vector. 
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FIG. 2: (color online) (a) The critical mode (the Laplacian 
eigenvector with ao = 422), (b) the activator pattern at the 
early evolution stage (t — 200), and (c) the stationary acti- 
vator pattern at the late stage (t = 1500). Nodes are ordered 
according to their degrees; with (d) showing the dependence 
of the degree on the node index. 

To perform numerical investigations of nonlinear dy- 
namics, a specific activator-inhibitor model is needed. 
For this purpose, we have chosen the Mimura-Murray 
model [ll|, which has been introduced to explain spa- 
tial nonuniformity of prey-predator populations in eco- 
logical systems. The chosen model is given by the 
equations f(U, V) = {(35 + 16U - U 2 )/9 - V}U and 
g(U,V) = {U - (1 + 2V/5)}V, whose fixed point is 
(U^°\V^) = (5,10). As an example of random net- 
works, we take below the Barabasi- Albert scale-free net- 
works [13[ of size TV = 1000. They are generated by the 
preferential attachment algorithm, starting from 10 fully 
connected initial nodes and adding m = 10 new connec- 
tions at each iteration step, so that their mean degree is 
(k) ~2m = 20. 

Figure [1] displays the real part of the linear growth rate 
Re A^ for three values of a with e = 0.12. The growth 
rate can become positive for a > a c = 15.5. The largest 
growth rate is reached at ao = 422 for a = 15.6. Near 
Im A^ = holds. In Fig. O we compare the 
critical mode (ao = 422) with two snapshots of the ac- 
tual activator patterns at the early and late stages of the 
evolution, as yielded by numerical integration of Eqs. (pQ) 
with e = 0.12 and a = 15.6. Here and below, node in- 
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dices {i} are sorted in the decreasing order of their node 
degrees {ki} so that k\ > &2 > • • • > kjy holds, which 
is useful in visualizing the Turing patterns exhibited by 
our system. 

Starting from almost uniform initial conditions with 
tiny perturbations, exponential growth is observed at 
the early stage. The activator pattern at this stage, 
Fig. [2(b), is similar to the critical mode, Fig. 0(a), with 
the deviations due to contributions from a few neighbor- 
ing modes that are already excited to some extent. Later 
on, however, strong nonlinear effects develop, and the fi- 
nal stationary pattern, Fig. [2(c), becomes very different 
from the one determined by the critical mode. Observ- 
ing the nonlinear development (see Video 0), we notice 
that some elements get progressively kicked off the main 
group near the destabilized uniform solution in this pro- 
cess. Eventually, in the asymptotic stationary state, all 
elements become separated into two groups. The corre- 
spondence between the indices and the degrees for the 
considered scale- free network is displayed in Fig. [2(d). 
The separation into two groups occurs only for the ele- 
ments on the nodes with relatively small degrees (roughly 
i > 200, hi < 24), while the elements on the nodes with 
high degrees (i < 200, ki > 24) do not separate. 
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FIG. 3: (color online) (a) Amplitude A of the Turing pat- 
tern vs. the diffusion ratio a; variation direction of a is indi- 
cated by arrows, (b) Stationary Turing patterns at parameter 
points P (a = 17.0), Q (a = 13.5), and R (a = 12.8). The 
inset in (a) shows the blowup near R. 

The outcome of nonlinear evolution depends sen- 
sitively on initial conditions. Different Turing pat- 
terns are possible at the same parameter values and 
strong hysteresis effects have been observed. As 
an example, Fig. [3(a) shows how the amplitude 
of the stationary Turing pattern, defined as A = 
E*=i {(Ui - U^) 2 + (Vi - V (0) ) 2 }, varies under grad- 
ual variation of the parameter a in the upward or down- 
ward directions. Stationary patterns observed at points 
P, Q, and R in Fig. [3(a) are presented in Fig. [3(b). As 
a was increased starting from the uniform initial condi- 
tion, the Turing instability took place at a = <r c , with the 
amplitude A suddenly jumping up to a high value that 
corresponds to the appearance of a kicked-off group. If 
a was further increased, the amplitude A grew. Starting 



to decrease a, we did not however observe a drop down 
at a = a c . Instead, a punctuated decrease in the ampli- 
tude A, which is characterized by many relatively small 
steps, was found. Reversing the direction of change of 
the parameter a at different points, many coexisting so- 
lution branches could be identified. The characteristics of 
Turing patterns vary with their amplitudes. When A is 
close to zero (point R in Fig. [3]), only a few kicked-off el- 
ements remain in the system. Because of the hysteresis, 
solutions with only a small number of destabilized ele- 
ments can coexist with the linearly stable uniform state 
in this region, so that the localized Turing patterns are 
found for a < <j c . 

The properties of the developed Turing patterns above 
the instability boundary (<r > <j c ) can be relatively well 
understood by using the mean-field approximation, pre- 
viously used for the analysis of epidemics spreadin g on 
networks [l4| and for networks of phase oscillators [15[. 
We start by writing Eqs. (pQ) in the form 



U i (t)=f(U i ,V i )+e(hP-k i U i ), 
V i {t)=g{U i ,V i )+ae{h { p -hVi), 



(4) 



where 



local 



fields felt 

(V) 



by each 



element, h 



AijVj, are introduced. 



These local fields are further approximated as h\ U ^ ~ 
kiHM and h { P ~ hH^, where global mean fields are 
defined by = Y,f=i ™jUj and = Y% =1 WjVj. 

The weights wj = kj/ (SjLi k 'j) = kj/k to tai take into 
account the difference in the contributions of different 
nodes to the global mean field, depending on their de- 
grees (cf. 0,[I3). Thus, the local fields are taken to be 
proportional to the degree of a node, ignoring the details 
of its actual connections. 

With this approximation, each element interacts only 
with the global mean fields, and its dynamics is described 

by 



U(t) = f(U,V) 
V(t)=g(U,V) 



U), 
-V). 



(5) 



We have dropped here the index z, since all elements obey 
the same dynamics, and introduced the parameter (3 = 
eki. If diffusion ratio a is fixed and the global mean fields 
H( u ) and are given, this parameter (3 plays the role 
of a bifurcation parameter that controls the dynamics of 
each element. Equations ([5|) have a single stable fixed 
point when (3 = (i.e. e = 0), and, as (3 is increased, this 
system undergoes imperfect pitchfork bifurcations that 
give rise to two new stable fixed points. 

We have computed stationary Turing patterns by nu- 
merical integration of Eqs. ([Tj) and determined the respec- 
tive global mean fields and at a = 15.6 and 
a = 30. Substituting these computed global mean fields 
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FIG. 4: (color online) Stationary Turing patterns compared 
with the bifurcation diagrams of a single element coupled to 
global mean fields. The parameters are e = 0.12 and (a,b) 
a = 15.6, (c,d) a = 30. Black curves (dots) indicate sta- 
ble branches, gray curves (dots) correspond to the unsta- 
ble branch. Crosses show the computed Turing patterns. 
The global mean fields are (H (U \H (V) ) = (4.95,9.97) for 
a = 15.6, and (if (t7) , if (v) ) = (4.8, 9.9) for a = 30. 



into Eqs. ([5]), bifurcation diagrams for a single Mimura- 
Murray element have been obtained (solid curves in Fig.[H 
(a,c)). These diagrams can be compared with the actual 
stationary Turing patterns. Each node i in the network is 
characterized by its degree fc^, so that it possesses a cer- 
tain value of the bifurcation parameter, (3 = eki. There- 
fore, the Turing pattern can be projected onto these bi- 
furcation diagrams, as shown by crosses in Fig. [H(a,c). 
We see a relatively good agreement between the stable 
branches and the data from the actual Turing patterns. 
Furthermore, we directly compare in Fig.[H(b,d) the com- 
puted Turing patterns with the mean- field predictions, 
based on Eqs. ([5]). We see that the Turing patterns are 
nicely fitted by the stable branches, though the scatter- 
ing of numerical data gets enhanced near the branching 
points. In this way, the fully developed Turing patterns 
in our system are essentially explained by the bifurcation 
diagrams of a single element coupled to constant global 
mean fields, with the coupling strength determined by 
the degree of the respective network node. 

Thus, we have constructed a general linear theory of 
Turing patterns in the activator-inhibitor systems on net- 
works and have also performed numerical investigations 
of nonlinear pattern formation. The Turing instability 
takes place on any networks as the ratio of diffusion con- 
stants of the inhibitor and activator species is increased. 
The critical mode corresponds to a certain eigenvector 
of the Laplacian matrix, but the final Turing patterns 
may be very different from the critical mode; multista- 
bility of solutions and hysteresis are effects are typically 



observed. Localized Turing patterns below the linear in- 
stability threshold have been found. 

To understand this behavior, it should be taken into 
account that, while being relatively large (N = 1000), 
the considered random networks had nonetheless only 
the small diameters of 4 and, moreover, each node in a 
network had 20 neighbors on the average. Therefore, dif- 
fusional mixing should be very strong in such activator- 
inhibitor systems, and behavior characteristic of small 
well-mixed spatial volumes can be expected; each ele- 
ment in a network may be viewed as interacting with 
some global mean fields. Indeed, separation of elements 
into two clusters has previously been reported for globally 
coupled activator-inhibitor systems under similar condi- 
tions [12]. The coupling to the global mean fields is 
however essentially heterogeneous in our case, because 
it is determined by the degrees of the respective network 
nodes which vary greatly. Although not included here, in- 
vestigations have also been undertaken by us for the clas- 
sical Brusselator model and for the Erdos-Renyi random 
networks, and similar results have been obtained. The 
constructed mean-field theory of nonlinear Turing pat- 
tern formation is applicable for various large random net- 
works of small diameters and different activator-inhibitor 
systems. 
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